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The single-particle density is the most basic quantity that can be calculated from a 
given many-body wave function. It provides the probability to find a particle at a given 
position when the average over many realizations of an experiment is taken. However, 
the outcome of single experimental shots of ultracold atom experiments is determined 
by the A/^-particle probability density. This difference can lead to surprising results. 
For example, independent Bose-Einstein condensates (BECs) with definite particle 
numbers form interference fringes even though no fringes would be expected based 
on the single-particle density [1—4]. By drawing random deviates from the iV-particle 
probability density single experimental shots can be simulated from first principles 
[1, 3, 5]. However, obtaining expressions for the A/^-particle probability density of 
realistic time-dependent many-body systems has so far been elusive. Here, we show 
how single experimental shots of general ultracold bosonic systems can be simulated 
based on numerical solutions of the many-body Schrodinger equation. We show how full 
counting distributions of observables involving any number of particles can be obtained 
and how correlation functions of any order can be evaluated. As examples we show 
the appearance of interference fringes in interacting independent BECs, fluctuations 
in the collisions of strongly attractive BECs, the appearance of randomly fluctuating 
vortices in rotating systems and the center of mass fluctuations of attractive BECs 
in a harmonic trap. The method described is broadly applicable to bosonic many- 
body systems whose phenomenology is driven by information beyond what is typically 
available in low-order correlation functions. 


Let us briefly outline how single experimental shots 
can be simulated from a general many-body wave func¬ 
tion At. The probability to find N particles at posi¬ 
tions ri,..., ftv in a many-body system is determined by 
the A-particle probability distribution P(ri,..., fat) = 
|^(ri,..., riv)P- In experiments on ultracold bosons 
snapshots of the positions of the particles are taken and 
single experimental shots sample the A-particle proba¬ 
bility distribution. This distribution is high-dimensional 
and sampling it directly from a given A-boson wave func¬ 
tion is hopeless. However, it can be rewritten as a prod¬ 
uct of conditional probabilities 

P(ri,..., Fat) = P(fi)P(f 2 |fi) x• • •xP(FAr|FAr-i,..., fi), 

( 1 ) 

where e.g. P(f 2 |fi) denotes the conditional probability 
to find a particle at F 2 given that another particle is at 
Fi. By drawing fi from P(f), F 2 from P(f|fi), fs from 
P(f|f 2 ,fi), etc., one random deviate of P(fi, ..., fat) 
is generated. Obtaining the conditional probabilities 
in (1) is a formidable combinatorial problem though, 
even for special cases [1, 5]. Here, we provide a gen¬ 
eral algorithm to simulate single shots from any given 
A-boson wave function |4^) = where |n) = 

|ni,..., um) are configurations constructed by distribut¬ 
ing N bosons over M orbitals <fi. We apply this algo¬ 
rithm to many-body states obtained by solving the time- 
dependent many-body Schrodinger equation numerically 
using the multiconfigurational time-dependent Hartree 
for bosons method (MCTDHB) [6-8]. This combination 
of many-body Schrodinger dynamics and sampling of the 
A-particle probability allows us to simulate single exper¬ 
imental shots from first principles in realistic settings, see 


Methods for the algorithm and details. 
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Figure 1. Interference of independent interacting condensates. 
Two independent, repulsively interacting condensates collide 
in an elongated trap. Shown is the single-particle density 
(left column) and random deviates of the iV-particle density 
(right column) at different times. In the overlap region in¬ 
terference fringes show up in the A-particle density, but not 
in the single-particle density. The results are obtained by 
solving the many-body Schrodinger equation in two spatial 
dimensions. Parameter values: N = 10000 bosons. Interac¬ 
tion strength A = 4.95. See text for details. All quantities 
shown are dimensionless. 

It is instructive to briefly review Bose-Einstein conden¬ 
sation. A many-boson state is condensed if its reduced 
single-particle density matrix has exactly one nonzero 
eigenvalue pi of order N [9]. The eigenvalues pi are 
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known as natural occupations, the eigenvectors as nat¬ 
ural orbitals. The BEC is fragmented if more than one 
eigenvalue pi is of order N [10, 11], see Methods for de¬ 
tails. Fully condensed states, i.e. states with pi = N^ 
are of the form (j)(Yi)(j)(Y 2 ) x • • • x 0 (rAr). ( 1 ) then be¬ 
comes a trivial product of independent, identical proba¬ 
bility distributions, and there are no correlations between 
particles. For instance, Gross-Pitaevskii (GP) mean-field 
states are of this form. Any other state, in particular 
fragmented states, exhibit correlations and vice versa any 
correlated state is to some degree fragmented. We will 
now show how fragmented BECs lead to macroscopically 
fluctuating outcomes in single shots. 

In the following we use dimensionless units h = m = 
1 and solve the time-dependent many-body Schrodinger 
equation using the MCTDHB method 

[ 6 - 8 ]. Here, 

N ^2 

H = ^---^ + V{ri) + XoY,Se{ri-rj) (2) 

i=l * i<j 

denotes a general many-body Hamiltonian in D dimen¬ 
sions with an external potential F(r) and a regularized 
contact interaction (^^(r) = (27re^)“^/^e“^ . We pa¬ 

rameterize the interaction strength by the mean-field pa¬ 
rameter A = Xo{N — 1), see Methods for details. 

Let us begin with an example of two interfering, in¬ 
dependent condensates of A" = 10000 bosons in an elon¬ 
gated trap with tight harmonic confinement along the z 
direction such that we can work in D = 2 dimensions and 
r = {x,y). We use F(r) = Vx{x) -\r Vy{y) + Vg{x) as an 
external potential, where Vx{x) and Vy{y) are harmonic 
traps and Vg{x) is an additional Gaussian potential that 
flattens the bottom of the trap along the x-direction. As 
an initial state we use two independent condensates, each 
of which is the mean-field ground state (corresponding to 
M = 1 in the MCTDHB formalism) of N/2 bosons of the 
displaced traps V±(r) = V{x^d,y) with d = 18.6 in har¬ 
monic oscillator units of the ^-direction at an interaction 
strength A = 4.95. The initial state |T(0)) = | A/ 2 , A/ 2 ) 
is fragmented with pi = p 2 = A/ 2 . We then solve 
the time-dependent many-body Schrodinger equation for 
|T(0)) using M = 2 orbitals. Fig. 1 shows the single¬ 
particle density as well as random deviates of the A- 
particle density at different times. The two condensates 
accelerate towards each other, collide and separate again. 
During the collision interference fringes appear in devi¬ 
ates of the A-particle density at locations that fluctuate 
randomly from shot to shot, but not in the single-particle 
density. This is also expected based on simplified models 
[1, 3]. However, here this result follows directly from the 
solution of the many-body Schrodinger equation. The 
interparticle interaction is weak here; interaction effects 
only become visible as ripples in the density after the 
collision and the natural occupations remain practically 
constant all along. 

We now go one step further and investigate collisions 
between strongly attractive independent condensates in 



Figure 2. Collision of independent attractively interacting 
condensates. Two independent attractively interacting con¬ 
densates collide in an elongated trap in two spatial dimen¬ 
sions. (a) Single-particle density at different times. The con¬ 
densates approach each other without spreading and bounce 
off one another, (b) Random deviates of the A-particle den¬ 
sity at the time of the collision. Correlations lead to either 
a single strongly localized density maximum containing prac¬ 
tically all particles or two smaller maxima containing about 
half the particles each, (c) Fragmentation of the condensate 
as a function of time. The initial state is two fold fragmented 
with pi/N = pn/N = 49.4%. During the collision two addi¬ 
tional natural occupations become significantly occupied and 
the system can no longer be separated into two independent 
condensates. Parameter values: A = 100 bosons. Interaction 
strength A = —5.94. See text for details. All quantities shown 
are dimensionless. 


the same trap. For this purpose we use A = 100 bosons 
at an interaction strength A = —5.94 which is about 2 % 
above the threshold for collapse of the GP mean-field 
ground state in this trap. For the initial state we first 
compute the many-body ground state of fifty bosons us¬ 
ing two orbitals and imaginary time-propagation. This 
ground state is highly condensed, pi/N = 98.7%. The 
initial state is then taken as the symmetrized product of 
the ground state and a displaced copy of it located at 
r = (—d, 0). Thus, the initial state has natural occupa¬ 
tions pi/N = P 2 /N = 49.4% and ps/N = P 4 /A = 0.6%. 
We then propagate this initial state using M = 4 orbitals. 

Fig. 2 (a) shows the single particle-density at differ¬ 
ent times. The condensates approach each other without 
spreading significantly, collide and separate again. Dur¬ 
ing the collision the single-particle density exhibits two 
maxima, the condensates seem to bounce off each other. 
However, single shots at the time of the collision reveal a 
different result, see Fig. 2 (b). In about half of all shots 
a strongly localized density maximum is visible, whereas 
in the other half two smaller well separated maxima ap¬ 
pear. We stress that at no point any type of (possibly 
random) phase relationship between the colliding parts 
was assumed. In fact, for independent condensates the 
assumption of a preexisting, but random relative phase 









3 



Figure 3. Fluctuating vortices. A repulsive condensate in 
the ground state of a harmonic trap is stirred by a rotating 
potential in two spatial dimensions. Over the course of time 
the system fragments and vortices appear at random posi¬ 
tions in single shots, (a) First column: single-particle density 
at different times. Second to fourth column: single shots at 
the same times, (b) Fragmentation of the condensate as a 
function of time. Starting from a condensed state, the sys¬ 
tem of bosons fragments as it is stirred. While the system 
is condensed single shots and the single-particle density look 
alike. When the system is fragmented vortices appear at ran¬ 
dom positions. Parameter values: N = 10000. Interaction 
strength: A = 17. See text for details. All quantities shown 
are dimensionless. 


is at variance with quantum mechanics [12]. The macro¬ 
scopic fluctuations in the outcomes follow directly from 
the intrinsic correlations of the many-body state. Fig. 2 
(c) shows the natural occupations of the system. As long 
as the condensates are far apart, the natural occupations 
remain close to their initial values. However, during the 
collision two additional natural orbitals become occupied 
indicating a buildup of even stronger correlations. As a 
consequence after the collision the system can no longer 
be separated into two independent condensates. 

In the previous two examples already the initial states 
were fragmented. We now turn to a system where frag¬ 
mentation builds up dynamically. Stirring a BEC can 
lead to fragmentation and vortex nucleation that cannot 
be explained within the mean-field framework of quan¬ 
tized vortices [5, 13]. Consider the ground state of a 
repulsively interacting BEC of N = 10000 bosons in a 
pancake shaped trap with = 1 at an interaction 

strength A = 17. We compute the many-body ground 
state using M = 2 orbitals which is practically fully con¬ 
densed with pi/N = 99.98%. We then switch on a time- 
dependent stirring potential Eg (r, t) = 
that imparts angular momentum onto the BEC. Here x{t) 
and y{t) vary harmonically and the amplitude r]{t) is lin¬ 
early ramped up from zero to a finite value and back 
down, see Methods for details. Eig. 3 (a) shows the den¬ 
sity together with single shots at different times. The 
evolution of the natural occupations is shown in Eig. 3 


Figure 4. Full counting distribution of the center of mass op¬ 
erator. Shown are 10000 random deviates of the center of 
mass operator of the ground state of an attractively inter¬ 
acting condensate in one spatial dimension. The center of 
mass fluctuations of the mean-field result (blue) are signifi¬ 
cantly smaller than those of the many-body results where the 
bosons are allowed to occupy M = 2, 3,10 (green, magenta, 
red) orbitals. The M = 10 result coincides with the exact 
analytical one (black). Parameter values: iV = 10 bosons; 
interaction strength A = —0.423, trap frequency cjx = 1/100. 
All quantities shown are dimensionless. 


(b). While the system is condensed, single shots repro¬ 
duce the single-particle density. Over the course of time 
an additional natural orbital becomes occupied and the 
BEC becomes correlated. As correlations build up the 
outcome of single shots fluctuates more and more and 
vortices appear at random locations in every single shot. 
This is in stark contrast to mean-field theory, where due 
to the lack of correlations vortices always appear at the 
same location. 

As a last example let us show how full distribution 
functions of A/'-body operators can be evaluated by sim¬ 
ulating single shots. Consider the ground state of N at¬ 
tractively interacting bosons in a harmonic trap, uJx = 
1/100, in one dimension, i.e. D = 1 and r = x. The 
exact wave function of the center of mass coordinate 
X = ^ J]) • Xi of the many-body ground state is given 

by a Gaussian ^^^(X) = with 

Xm}) = 1 / \/Nujx [14]. On the other hand, the mean- 
field ground state is uncorrelated and hence its center 
of mass width is given by X^f = where 

= ( 0 m/ 1 ^^ 10 m/) is the variance of the mean-field 
orbital 0m/, see Methods. In the limit of a weak trap, 
ujx ^ the mean-field solution approaches a soliton 
with (Tmf = 7 r/(v^|A|). Thus, for sufficiently strong at¬ 
tractive interaction exceeds X^j. We compute the 
ground state of X = 10 bosons at an interaction strength 
A = —0.423 using imaginary time-propagation for differ¬ 
ent numbers of orbitals. Erom the obtained ground states 
we generate 10000 random deviates of the center of mass 
coordinate. Eig. 4 shows fits to the obtained histograms 
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of the center of mass deviates together with the exact 
center of mass distribution. The many-body result for 
M = 10 orbitals is indistinguishable from the exact one 
and significantly broader than the mean-field (M = 1 ) 
result. In the present example the many-body correla¬ 
tions are the cause for the onset of the delocalization of 
the ground state. 


METHODS 


Bose-Einstein condensation. 

For an A^-boson state |T) = Cfi{t)\n) and a bosonic 
field operator T(r) = ^jbj(f)j{r) the reduced single¬ 
particle density matrix is defined as 

p(i)(r|r') = (4'|4't(r')^(r)|«') = (3) 

with pij = By diagonalizing pij one ob¬ 
tains p^^^(r|r') = The eigenval¬ 

ues Pi > P 2 > • • • are known as natural occupations, the 
eigenvectors as natural orbitals. If there is only 

one eigenvalue pi = 0{N) the BEG is condensed [9], if 
more than one pi = 0{N) the BEG is fragmented [10, 11]. 
The diagonal p(r) = p^^^(r|r' = r) is the single-particle 
density of the N-hoson wave function. 


Single Shot Algorithm. 

Here we show how single shots can be simulated from 
a general A-boson wave function expanded in M or¬ 
bitals |T) = where \n) = |ni,...,nM) and 

^ 2=1 Special cases (for M = 2 ) have been 

treated in earlier works [1, 5]. The goal is to draw the 
positions ri,..., r^v of A bosons from the probability dis¬ 
tribution P(ri,..., Fat). We achieve this by evaluating 
the conditional probabilities in (1). Eor this purpose we 
define reduced wave functions 

|^(/c)) ^ iik-0 

' ^ \A4T(r,)|T(^-i)), if A: = 1,...,7V-1 

of n = N — k bosons with normalization constants 
A4- The respective single-particle densities are given by 
pk{r) = and A4 = /9fe-i(rfe)“i/^. 

The first position ri is drawn from P(r) = po{r)/N. 
Assuming that positions r/.,..., ri have already been 
drawn, the conditional probability density for the next 
particle P(r|r/e,..., ri) = P(r, r/^,..., ri)/P(r/c,..., ri) 
is given by 


P(r|rfc,...,ri) cx/9fe(r), (5) 

since P(r/c,..., ri) is a constant. The problem is 
thus reduced to obtaining the wave function = 


from the wave function = 

C? where the sums over run over all config¬ 

urations of n and n 1 bosons, respectively. Defining 
= (ni,..., rig -h 1,..., um) one finds from (4) 

M 

CP = A4 y] </>g(r)ciy xAT+I (6) 

9=1 

Using ( 6 ) in a general M orbital algorithm requires an 
ordering of the configurations |fi) for all particle 

numbers n = 1 ,..., A. Gombinadics [7] provide such an 
ordering by associating the index 




M-l 

^m) = 1 + 

2=1 


n + M ■ 


1-i- 
M-i 





( 7 ) 

with each configuration In). Using (7) all coefficients 
can then be obtained by evaluating the sums in ( 6 ) and 
A4 is determined by normalization. Using the coefficients 
we evaluate p/c(r) and by means of (5) we then draw 
r/e+i from P(r|r/c,..., ri). This concludes the algorithm 
to simulate single shots. It is now easy to see that also 
correlation functions of arbitrary order can be evaluated. 
By realizing that 


k 

i=i 

(8) 

the k-th order correlation function is evaluated at 
ri,..., r/c as the product of the reduced densi¬ 
ties pj-i{Tj). To evaluate the correlation function 
(T|T’f'(ri)... T’f'(r/c)T(r/c)... T(ri)|T) the only modifi¬ 
cation to the single shot algorithm above consists in 
choosing the positions ri,..., r/. rather than drawing 
them randomly. 


MCTDHB. 

In the MGTDHB [ 6 - 8 ] method the many-boson wave 
function is expanded in all configurations that can be 
constructed by distributing N bosons over M time- 
dependent orbitals 0i(r, t). The ansatz for the time- 
dependent many-boson wave function reads: 

= '^Cfi{t)\n;t) (9) 

n 

In (9) the Cfi{t) are time-dependent expansion coeffi¬ 
cients and the jn; t) are time-dependent permanents built 
from the orbitals 0i(r, t). The MGTDHB equations of 
motion are derived by requiring stationarity of the many- 
body Schrodinger action functional 

- Ylk,j=i - 4i]}, (10) 
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with respect to variations of the coefficients and the or¬ 
bitals. The jJikj{t) are time-dependent Lagrange multi¬ 
pliers that ensure the orthonormality of the orbitals. For 
bosons interacting via a delta-function interaction and 
M = 1 the MCTDHB equations of motion reduce to the 
time-dependent Gross-Pitaevskii equation. For more in¬ 
formation see the literature [6-8]. 


Parameters. 

For the D = 2 dimensional simulations in this work 
we assume tight harmonic confinement with a frequency 
ujz and a harmonic oscillator length Iz = yh/{rruJ^ 
along the z -direction. The bosons interact via a reg¬ 
ularized contact interaction potential with 

^e(r) = and a dimensionless interac¬ 

tion strength Aq = V^a/lz^ where a is the scattering 
length and m the mass of boson. We note that it is im¬ 
portant to regularize contact interaction potentials for 
D 1 [15^ 16]. J-'he contributions to the external po¬ 
tential are given by Vx{x) = 

and Vg{x) = Ce~^ , with C = We obtain 

dimensionless units h = m = 1 and the Hamiltonian (2) 
by measuring energy in units of hujy^ length in units of 
ly = ^Jh/{mujy) and time in units of 1/ujy. We use a 
plane wave discrete variable representation to represent 
all orbitals and operators. The width of the contact inter¬ 
action is e = 0.15 and the grid spacing is Ax = Ay = e/2 
for all simulations in this work. For the elongated trap 
the parameter values are = 0.07, = 1 and a = 10 


on a grid [—43.2,43.2] x [—3.6,3.6]. For the rotating 
BEG the parameter values are = ^y = ^ and r]{t) 
is linearly ramped up from zero to r]max = 0.1 over 
a time span tr = 80. r]{t) is then kept constant for 

tup = 220 and ramped back down to zero over a time 
span tr. The potential W(r, t) = ^r]{t)[x{t)‘^ — y{t)‘^] ro¬ 
tates harmonically with x{t) = xcos(flt) + ysm{Qt) and 
y{t) = —xsm{ftt) -h^cos(flt), where = 7 r/ 4 . The grid 
size is [— 8 , 8 ] x [— 8 , 8 ]. 

For the D = 1 dimensional simulations we as¬ 

sume tight harmonic confinement along the y- and 
^-directions with a radial frequency uj±_ = ooy = ooz 
and an oscillator length l±_ = y^h/{muj±). The contact 
interaction potential is then given by with 

Se{x) = We use as the unit of 

energy and l± as the unit length. The dimensionless 
interaction strength is then given by Aq = 2a/l_\_. The 
harmonic potential along the x-direction = 1/100 
is much weaker than the radial confinement uj± = 1 . 
The grid size is [—90,90]. The Gross-Pitaevskii soli- 
ton solutio n on an infinite line takes on the form 
4>mf{x) = yA/4 sech(Aa;/2). 

A. Image processing. 

The histograms of the positions of particles obtained 
using the single shot algorithm have a resolution that is 
determined by the grid spacing. For better visibility and 
in analogy to a realistic imaging system we convoluted 
the data points of each histogram with a point-spread 
function (PSF). As a PSF we used a Gaussian of width 
3x3 pixels. 
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